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Interpolating  Polynomial  Macro- Elements 
with  Tension  Properties 


Paolo  Costantini  and  Carla  Manni 


Abstract.  In  this  paper  we  present  the  construction  of  nine-parameter 
polynomial  macro  elements,  based  on  the  classical  Powell-Sabin  split, 
which  can  be  connected  to  form  a C1  surface.  Variable  degrees,  which 
act  as  independent  tension  parameters,  are  associated  with  any  vertex  of 
the  triangulation,  i.e.  to  any  interpolation  point. 


§1.  Introduction 

Among  the  several  approaches  for  avoiding  extraneous  inflection  points  in 
interpolating  functions,  the  so  called  tension  methods  are  the  oldest  and  prob- 
ably the  most  famous  ones.  Basically  they  consist  of  Ck,  k > 1,  piecewise 
functions  depending  on  a set  of  parameters  which  are  selected  in  a local  or 
global  way  to  control  the  shape  of  the  interpolants,  stretching  their  patches 
between  data  points. 

Piecewise  polynomial  splines  with  variable  degrees  have  turned  out  to  be  a 
useful  alternative  to  classical  (exponential  or  rational)  tension  methods,  and 
have  been  successfully  applied  both  in  free-form  design  and  in  constrained 
interpolation  of  spatial  data  [1],  On  the  other  hand,  a limited  number  of 
methods  for  constrained  interpolation  of  bivariate  scattered  data  are  at  present 
available,  and  very  few  of  them  offer  the  possibility  of  controlling  the  shape 
via  tension  parameters  (see  for  example  [2,3]  and  references  quoted  therein). 

Let  a set  of  scattered  points  ( Xi,yi,fi ),  i = 0,1, be  given,  and 
suppose  they  have  been  associated  with  a proper  triangulation  T . We  are 
interested  in  local  methods  where  the  polynomial  pieces  of  a spline  are  de- 
termined one  triangle  at  time  using  only  local  data.  Such  methods  are  called 
macro-element  methods.  The  aim  of  the  paper  is  to  describe  a new  class  of 
variable  degree  polynomial  triangular  macro-elements,  and  to  show  their  ten- 
sion properties.  These  polynomial  patches,  which  are  based  on  the  classical 
Powell-Sabin  split  of  a triangle,  can  be  connected  to  form  a C1  interpolating 
function. 
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Typically,  the  main  drawback  of  macro-element  methods  is  in  the  strong 
influence  of  the  triangulation  on  the  shape  of  the  interpolating  surface,  and, 
as  far  as  we  know,  it  is  still  not  clear  how  to  construct  a “good”  triangulation. 
Given  that  any  interpolating  scheme  based  on  triangular  macro-elements  can- 
not be  completely  independent  of  the  triangulation,  we  can  nevertheless  try 
to  reduce  this  dependence.  From  this  point  of  view,  the  main  advantage  over 
existing  methods  is  in  the  fact  that  we  have  a variable  degree  associated  to 
each  vertex  of  the  triangulation;  therefore  the  modification  in  the  shape  of 
the  interpolant  is  similar  for  all  the  patches  around  the  same  common  vertex 
(see  Fig.  8 right).  In  other  words,  although  constructed  over  T,  the  tension 
parameters  are  more  related  to  the  interpolation  points  than  to  the  edges  of 
the  adopted  triangulation. 

It  is  worthwhile  to  anticipate  that  the  possibly  high  degrees  we  use  for  our 
construction  act  only  as  tension  parameters,  and  do  not  modify  the  basic  struc- 
ture of  the  macro-element.  In  other  words,  we  always  use  a nine-parameter 
macro-element,  and  the  computational  complexity  is  almost  independent  of 
the  size  of  the  degrees  used. 

The  scheme,  being  local,  requires  that  the  gradients  are  also  known  at  the 
vertices  of  the  triangulation.  If  this  information  is  unavailable,  gradients  can 
be  recovered  from  the  data  points.  In  the  presented  numerical  test,  gradients 
have  been  computed  from  the  data  according  to  the  classical  least  square 
strategy. 

The  paper  is  divided  into  six  sections.  In  the  next  section  we  introduce 
some  notations.  Sections  3 and  4 are  devoted  to  the  construction  of  the  control 
points  defining  the  macro-element,  whose  properties  are  briefly  discussed  in 
Section  5.  We  end  with  Section  6 where  a graphical  example  is  presented. 


§2.  Notation  and  Preliminaries 

In  this  section  we  introduce  some  notation.  To  aid  in  comprehension,  we 
notice  that  points  and  vectors  in  1R2  and  in  1R3  have  been  denoted  by  bold- 
faced characters  unless  classical  notations  (as  for  gradients)  have  been  used. 
As  usual,  we  describe  the  polynomial  macro-element  in  terms  of  its  Bernstein- 
Bezier  control  points.  Let  Pr,  r = 1, 2, 3,  be  three  non-collinear  points  in  1R2, 
and  let  T denote  the  triangle  they  form.  An  n-degree  Bernstein  polynomial 
has  the  form 


b(x,y;n)  = b(u,v,w,n)  := 


E 

i+j+k-n 


7l\ 

i\j\k\ 


lijk  'll  11)  i 


where,  setting  P = [rr , y]T  € 1R2,  u — u(x,y),  v — v(x,y),  w = w(x,y),  are 
the  barycentric  coordinates  of  P with  respect  to  the  vertices  of  T,  that  is 

P = wPi  + vP2  + WP3;  u + v + w — 1, 

and  lijk,  i + j + k = n,  are  the  Bernstein  ordinates  of  b(., ., .;  n). 
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Fig.  1.  Powell-Sabin  split  of  a triangle. 


Setting  x = x(u,v,w),  y = y(u,v,w),  the  points  in  ]R3 
-x(i  i k) l 

' n ’ n ’ n ' 


Jijk 


v(-  1 -1 

y\n'  n'  n> 

hjk 


i,j,k  > 0 , i + j + k = n, 


are  called  control  points  [4]. 

We  finally  recall  the  so-called  Powell-Sabin  split  ([5])  of  a given  triangle, 
T,  which  consists  of  dividing  T :=  PiP2P3  into  six  mini-triangles  (Fig.  1): 

T<x’°>  :=  P2M!R,  T(2’0)  :=  P3M2R,  :=  P1M3R, 

T(1(1)  :=  M1P3R,  :=  M2PiR,  T P-1)  :=  M3P2R, 


where 

R = /3lPl  + 02^2  + P3P3  ) 01+0 2 + /?3  = 1> 
is  a point  internal  to  T and 


Mj  = (1  - ai) Pi+i  + a,P,+2,  0 < a,  < 1, 

is  a point  internal  to  the  edge  of  T opposite  to  Pj.  Here,  and  in  the  following, 
indices  will  be  considered  modulus  3. 

We  will  denote  by  L L^.9  the  control  points  of  Bernstein  polynomi- 
als in  the  mini  triangle  Tl'p'qi  (see  below).  Let  the  data 

(Pi,/i  = /(Pi),v/i  = v/(Pi)),PieR2,  * = 1,2,3,  (l) 

be  given.  As  mentioned  in  the  introduction,  our  goal  is  to  construct  a C 1 
polynomial  macro-element  on  T interpolating  the  data  and  having  tension 
properties. 

The  macro-element  will  be  obtained  considering  a Powell-Sabin  split  of 
T,  and  constructing  in  each  mini  triangle  T+’q>  a Bernstein  polynomial  via 
suitable  control  points,  L . To  obtain  the  final  control  points  (FCP)  , 
we  follow  two  basic  steps:  first  we  construct  the  basic  control  points  (BCP) 
L , then  we  modify  them  to  reach  the  required  smoothness  of  the  macro- 
element. 
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§3.  Defining  the  Basic  Control  Points 

For  each  vertex  P,,  let  us  consider  an  associated  given  degree 

7i;  >3,  71;  € IN. 


Let  us  now  describe  the  construction  of  the  BCP  considering  for  the  sake  of 
simplicity  only  the  mini  triangle  T'3,0^.  See  Fig.  2 for  the  role  of  indices. 

First  of  all  we  assume  interpolation  conditions  (for  the  position  and  the 
gradient)  at  Pi  (see  in  Fig.  3 left)  that  is 

,(3,o)  f 

‘"i.O.O  — J 1 > 

^-)i.i.o==/i  + ^^/i,PiM3),  (2) 

tt, o,i  :=/i  + ^-(V/1;PiR). 

’ Tl\ 

In  order  to  define  the  BCP  around  M3,  let  us  consider  the  univariate  piecewise 
linear  function,  l3,  defined  along  the  edge  P1P2  having  breakpoints  at 

Pi,  P1  + -P1M3,  p2--m3p2,  p2, 

Til  n2 

interpolating  / and  its  derivatives  at  the  extremes  of  the  edge.  We  define  (see 
“®”  in  Fig.  3 left) 

C.o  :=  h(M3), 

lt°Uo  ■■=  Is  (m3  - ^PiM3)  , (3) 

C-1.1  :=  *3(M3)  + ±d3, 

71 1 


where 


d3  :=  ((1  - a3)V/(Pi)  + a3V/(P2),  M3R) . 

Moreover,  we  require  that  Lj^_21  (see  “©”  in  Fig.  3 left)  belongs  to  the 
plane  through 


(3,0)  T (3,0)  T (3, 

l,ni-l,0i  ^O.nj.O’  -^O.ni-l.l 


Concerning  the  remaining  control  points,  we  assume  that  the  not  yet  defined 
control  points  of  the  first  two  rows  in  the  mini  triangle  parallel  to  PiP2  (see 
“o”  in  Fig.  3 left)  belong  to  the  straight  lines  through  the  above  defined  control 
points,  that  is 

t (3,o)  _i-l,  (3,0)  , ("1  - 2)  - (i  - 1)T  (3,0) 

Kni-z.O  ni  _ 2L‘n1-l,hO  + ni-2  i, «i-l,0  > 

T (3,0)  * - 1 T (3,0)  , ("1  - 2)  “ (*  “ 1)  T (3,0) 

-ni  _ 2^1-1, 0,1+  m — 2 ’ 

i =rii  — 2, ...  ,2  . 


(4) 
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Fig.  3.  The  construction  of  the  BCP  for  the  split  of  Figure  1 with  n\  = 5,  712  = 
3,  T13  = 8.  Left:  projection  of  L in  the  x,  y plane;  are  deter- 
mined via  interpolation  conditions  (2);  “®"  by  (3);  “o”  by  (4)  and  “o” 
by  coplanarity  conditions.  Right:  , p = 1,2  3,  q = 0, 1.  Note  the 

discontinuities  across  the  edges  M,R. 


Moreover,  we  require  that  conditions  for  C 1 continuity  across  the  edge  PjR 
hold  [4],  Then,  in  particular, 

T (2,1)  T (2,1)  t (2,1)  T (3,0)  T (3,0)  t (3,0) 

1Jl,n1-2,l’lj0,n1-l,lilj0,ni-2,2  ~ ^m-2,0,2!  un-i -1,0, 1>  ±Jm-2,l,l> 

lie  onto  the  same  plane. 
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Fig.  4.  Left:  the  control  points  L jj'*  after  the  first  step,  n = 8.  Right:  The 
final  control  points  L after  the  second  step.  L:7’’^  as  in  Fig.  3. 


Finally,  the  central  control  points  L^’^,  k > 2,  (see  “o”  in  Fig.  3 left) 
are  assumed  to  lie  onto  the  plane  through 


(3,0)  T (1,0)  T (2, 
'm-2,0,2!  2,0,2’  ljn3-2,0,2 


Similarly  we  define  the  BCP  in  the  other  mini  triangles. 

The  above  construction  provides  a polynomial  macro-element  which 
turns  out  to  be  of  class  C 1 across  the  interior  edges  P,R. 

Moreover,  we  assume  that  the  points  M lying  on  the  internal  edges  of  the 
initial  triangulation  T also  lie  on  the  straight  lines  joining  the  R points  of  the 
triangles  which  those  edges  separate  (see  Fig.  8).  This  classical  requirement, 
(2),  (3),  (4),  and  the  geometry  of  the  Powell-Sabin  split  ensure  C1  continuity 
of  two  macro-elements  across  the  boundary  edge  P,Pt+]  (see  Fig.  3 right). 

On  the  other  hand,  the  BCP  do  not  produce  in  general  a continuous 
macro-element  across  the  edges  M,R  unless  the  degrees  n}  are  equal  (see 
Fig.  3 right).  In  order  to  obtain  a C1  macro-element  without  imposing  any 
conditon  on  the  degrees,  we  modify  the  constructed  BCP.  This  will  be  de- 
scribed in  the  next  section. 


§4.  Obtaining  a C1  Macro  Element 

In  this  section  we  describe  how  to  modify  the  BCP  to  obtain  the  final  control 
points  (FCP)  producing  a C1  macro-element.  The  modified  FCP  will  be 
basically  obtained  from  the  BCP  via  the  degree-raising  process  in  two  steps. 
As  a first  step,  for  each  mini  triangle  T<-P:q'>  let  us  compute 

I . *.  3,  k > 0 , i+j  + k = n, 

the  control  points  obtained  from  L-7^,  i,j,  k > 0 , i + j + k = np+q+ 1,  by 
the  degree-raising  process  ([4])  from  the  degree  np+q+ 1 (that  is  the  degree 
associated  with  the  mini  triangle  T<-P'q'1)  to  the  degree  n :=  maxjni,  n2,  n3} 

(see  Fig.  4 left).  The  control  points  L allow  us  to  express  the  polynomial 
of  degree  np+q+i  defined  by  L,-7’’1^  as  a Bernstein  polynomial  of  degree  n. 
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We  emphasize  that,  due  to  the  geometry  of  the  split,  the  control  points 
define  a macro-element  which  is  of  class  C1  across  the  interior  edges 
MPR  if  and  only  if  the  control  points 

f-(p.o)  w(p> o)  _ w(p,i)  w(p,i)  . _ n ^ i 

^0 ,n-k,k  ~ "n-fc-l.l.fc)  * U,  • • ■ ,Tl  i, 

are  collinear.  Moreover,  we  notice  that,  due  to  the  construction  of  the  BCP 
around  Mp  and  near  R and  to  the  properties  of  the  degree-raising  process, 
the  control  points 

w(p> °)  f-(p.o)  _ j-Cp.1)  w(p,i)  , _ n i _ i 

■^1,71— 1— k,ki  -^0,71— fc,fc  k,0,kf  ^n  — k— 1,1, fe)  * U,  1,71  1, 


lie  on  the  same  plane,  then  they  are  collinear  due  to  the  geometry  of  the  split 
(see  Fig.  4 left).  Then  in  order  to  obtain  C 1 continuity  across  the  edge  MpR, 
we  simply  consider  a second  step  in  which  we  modify 


w(p.o)  w(p.l)  , _ 0 


,n  2,  p — 1, 2, 3, 


(5) 


imposing  that  (see  Fig.  4 right)  they  lie  on  the  segment  through 


T (P'0)  r (Pd)  r n 


n — 2,  p = 1,2,3. 


§5.  Properties  of  the  Macro— Element 

In  this  section  we  analyze  the  interpolation,  smoothness  and  tension  properties 
of  the  macro-element  defined  by  the  FCP  constructed  in  Section  4.  First  of  all, 
we  notice  that  the  construction  of  the  macro-element  is  completely  local:  it 
only  depends  on  the  data  at  the  vertices  of  T and  on  the  degrees  nj  associated 
with  the  vertices  which  are  given  input  parameters. 

Theorem  1.  The  polynomial  macro-element  defined  by  the  control  points 

I'ij.'fc  . P = !>  2>  3,  q = 0, 1,  i,j,k>  0 , i + j + k = n, 

interpolates  the  data  (1)  and  is  of  class  C1(T).  Moreover,  let  T be  a given 
triangulation  equipped  with  a classical  Powell-Sabin  split,  and  with  a fixed 
degree  associated  with  any  vertex.  Then  the  collection  of  the  macro-elements 
corresponding  to  the  triangles  ofT  produces  an  interpolating  surface  of  class 
C1. 

Proof:  The  BCP  defined  in  Section  3 produce  an  interpolating  macro- 
element which  is  of  class  C1  across  the  interior  edges  P;R  and  the  boundary 
edges  PjPj+i.  After  the  degree-raising  process,  only  the  control  points  in 
(5)  are  modified  in  order  to  obtain  C 1 continuity  across  M,R.  Since  the  C1 
continuity  across  one  edge  only  depends  on  the  control  points  lying  on  the 
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m = n2  = n3  = 3 


ni  = 3,  n2  = n3  = 6 


m = n2  = n3  = 15 


Fig.  5.  The  local  tension  effect  of  the  degrees:  vertices  are  numbered  counter- 
clockwise from  the  origin. 

first  two  rows  parallel  to  that  edge  ([4]),  the  second  step  in  Section  4 does  not 
affect  Cl  continuity  across  P,R  and  PiPi+i  □ 

As  mentioned  before,  the  degrees  nj  are  free  input  parameters.  From 
the  construction  described  in  Section  3,  it  is  clear  that  increasing  their  val- 
ues causes  the  BCP  to  approach  the  plane  interpolating  (P i = 1,2,3. 
Similarly,  the  FCP  approach  the  same  plane  because  they  have  been  obtained 
via  a degree-raising  process,  that  is  via  a convex  combination  of  the  BCP. 
Therefore,  the  same  property  is  shared  by  the  macro-element,  due  to  the  con- 
vex hull  property  of  the  Bernstein  representation.  We  summarize  the  tension 
properties  of  our  macro-element  with  the  following  theorem. 

Theorem  2.  If  n\ , ri2,  n$  — > Too,  then  the  polynomial  macro-element  de- 
fined by  the  control  points 

P = 1.2,3,  q = 0, 1,  i,j,  k>  0 , i + j + k = n, 

approaches  the  plane  through  (P,,  /,),  i = 1,2,3. 

We  end  this  section  by  emphasizing  that  the  degrees  act  as  local  tension 
parameters:  each  degree  affects  the  shape  of  the  interpolating  surface  only 
around  the  associated  vertex  (see  Fig.  8 right),  and,  as  previously  said,  this  is 
the  main  feature  of  this  method.  The  increase  of  a degree  pushes  the  surface 
to  the  piecewise  linear  interpolant  around  the  corresponding  point,  giving  it 
a cuspidal  appearance.  This  local  tension  effect  is  clearly  shown  in  Fig.  5. 
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§6.  Numerical  Results 

In  this  section  we  present  a classical  graphical  example  to  show  the  perfor- 
mances of  the  macro-element  and  the  local  tension  effect  of  the  degrees. 

The  data  (Fig.  6 left)  are  taken  from  [6]  (see  also  Fig.  8 left  for  the  consid- 
ered triangulation).  The  interpolating  surfaces  obtained  by  using  the  proposed 
macro-element  with  uniform  degrees  rij  = 3 is  depicted  in  Figure  6 right.  Fig- 
ure 7 left  shows  the  interpolating  surfaces  obtained  by  using  uniform  degrees 
n,j  = 9.  The  tension  effect  due  to  the  increased  values  of  the  degrees  is  evi- 
dent, but  obviously  uniformly  distributed  over  all  of  the  surface.  On  the  other 
hand,  the  surface  in  Figure  7 right  has  been  obtained  considering  all  the  de- 
grees equal  to  3,  except  those  associated  with  the  five  vertices  as  depicted  in 
Figure  8 left.  The  local  tension  effect  of  the  degrees  is  clear.  The  influence 
region  of  the  increased  degrees  can  be  also  seen  in  Figure  8 right. 
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